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Abstract 

We test here the first stage of a route of modifications to be applied to the public GADGET2 code for dynamically 
identifying accretion centers during the collision process of two adjacent and identical gas cores. Each colliding core 
has a uniform density profile and rigid body rotation; its mass and size have been chosen to represent the observed 
core L 1544; for the thermal and rotational energy ratios with respect to the potential energy, we assume the values 
a = 0.3 and /3 = 0.1, respectively. These values favor the gravitational collapse of the core. We here study cases 
of both -head-on and off-center collisions, in which the pre-collision velocity increases the initial sound speed of the 
barotropic gas by up to several times. In a simulation the accretion centers are formed by the highest density particles, 
so we here report their location and properties in order to realize the collision effects on the collapsing and colliding 
cores. In one of the models we observe a roughly spherical distribution of accretion centers located at the front wave 
of the collision. In a forthcoming publication we will apply the full modified GADGET code to study the collision of 
turbulent cores. 

Keywords: stars: formation- physical processes: gravitational collapse, hydrodynamics- methods:numerical. 


1 Introduction 

Stars form out of collapsing prestellar cores, located usually in larger clumps that could be isolated or part of a larger 
gas filament. Massive stars could also form as a result of the collapse of a relatively large core, | [Zinnecker & Yorke (2007)| . 
Besides, the massive star formation scenario by means of gas collisions has gained observational support from recent 
space observation projects like ALMA, NANTEN, and the Spitzer telescope, all of which provide evidence of the 
occurrence of collisions among gas clouds and/or clumps, see section 1.1 of [Takahira et al.(2014)j and references 
there in. 

Besides, [Testi et al.(2000)| have presented high-resolution interferometric and single-dish observations of molec¬ 
ular gas in the Serpens cluster-forming core. It appears that the Serpens core encompasses at least three substructures, 
the northwest and southeast sub-clusters and the near-IR cluster. Thus, the spatially distinct structures are also kine¬ 
matically separated, while there is reasonable internal velocity coherence. Star formation is currently occurring 
simultaneously in the northwest and southeast sub-clusters, both of which contain roughly equal fractions of pre¬ 
stellar, proto-stellar, and IR sources. These facts point to the idea that Serpens was formed by means of a collision 
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between two or three very similar gas clumps. |Graves et al. (2010)] and |Duarte-Cabra l et al, (2010 )] have recentely 
observed the Serpens molecular cloud by using HARP at the James Clerk Maxwell Telescope and reported evidence 
that Serpents consists of two colliding sub clouds. 

Furthermore, |Furukawa et al. (2009)| , |Torii et al.(2011)| and |Churchwell et al. (2006)| , have reported about 
more recent observational evidence of cluster formation triggered by cloud-cloud collisions. They included the cases 
of RCW49, Westerlund2 and NGC3603 as examples. We then assume that collisions between smaller gas structures 
are also possible. 

Turbulence makes a big change in the spatial distribution of the interstellar gas as turbulent clouds become fil¬ 
amentary and flocculent and their gas is locally compressed such that dense cores can form and eventually become 
gravitationally unstable to form stars, see for instance |Federrath & Klessen (2012)| and |Padoan et al.(2014)|. It is 
within highly turbulent molecular clouds where core collisions may occur because of the random nature of the gas 
velocity favors the formation of small over-densities across the cloud, which might be already collapsing or at least 
about to start collapsing. 

Besides, [Motte et al. (1998)) have shown evidence that several parts of the Oph cloud have been compressed by 
an external pressure, such as a supernova explosion, an expanding HII region, or stellar wind shells. [Scoville et al. (1986)| 
have also presented observational evidence of gas compression at the interface of two colliding molecular gas regions. 

Thus, there is sufficient observational evidence to support the idea that gas collisions in all their possible size struc¬ 
tures can occur. One can then ask what could be the influence of core collision on star formation. 

Collision between relatively smaller clumps of gas and larger clouds have been studied numerically for more 
than three decades, for instance, |Hausman (1978)) and {Lattanzio et al. (1985)] ; particularly, head-on collisions 
were considered by |Kimura & Tosa (1996)| , |Klein & Woods (1998)| and |Marinho & Lepine (2000)| . Some of 
the earliest works were performed with only a few thousand particles and therefore suffered from a lack of reso¬ 
lution. However, in the recent past, simulations have been developed with much higher spatial resolution, for instance 
|Burkert&~ AI ves (2009)| and (Anathpindika (2009a) |. The case of colliding gas structures starting from hydrody- 
namical equilibrium was considered by |Kitsionas & Whitworth (2007) | and | Anathpindika (2009b)|; [Anathpindika (2010) | 
has considered collisions between dissimilar gas structures. |Gomez et al. (2007)| and |Vazquez-Semanedi et al. (2 007) | 
and other authors as well, have studied the generation of turbulence and star formation at the shock front of head-on 
collisions. 

We here present numerical high resolution three dimensional hydrodynamical simulations of two rotating gas core 
collisions, including the self-gravity of the gas. We calculate all collision models with a modified GADGET2 code 
(which is describe in Section |3.3fr in order to describe the formation of accretion centers in a collision system. We 
must clarify that accretion centers are conceptually simlar but not identical to the sink particles. 

The sink particle technique, as implemented by |Ba te et al. (1995)| and jFederrath et al. (2010) |, has been very 
useful, as one of the most important issues in star formation numerical simulations is that of identifying peak density 
gas parcels and monitoring their evolution as closely as possible. In the case of the particle-based codes used to follow 
the gravitational collapse of cores, all the higher density particles are very closely packed, so that this set of particles 
is likely to be identified with a protostar. The sink particles play the role of proto-stars in the simulations and their 
physical properties are easily investigated in this scheme, which hopefully could be compared with observations. For 
instance, (Clark & Bonnell (2006)| have used the sink particle technique in their study of the mass function of cores 
by means of colliding clumpy flows. 

We devote this paper to making new numerical simulations of the collision process of two rotating cores, taking 
advantage of both the new capabilities of modern computers and improved numerical techniques. We emphasize that a 
core very similar to the one which we take here, was originally suggested by |Whithworth & Ward-Thompson (2001)| 
as an empirical model to study the collapse of the well-observed core L1544. The physical properties of this core fit 
more or less the observed average properties of the cores reported by (Tafalla et al.(2004)| and (Jijina et al.(1999)| . 
[Tafalla et al.(2004)) have presented a multi-line and continuum study of LI544, and according to them this core 
presents evidence to be in an advanced state of evolution towards the formation of one or two low mass stars. 

In this paper we take only the core size and mass as given by |Whithworth & Ward-Thompson (2001)| , and we 
additionally consider the core to be in counterclockwise rigid body rotation around the Z axis. We also include a 
density perturbation on the initial core model in order for the isolated collapsing core to develop a filament which 
eventually would fragment into a binary protostar system. 

We here study collision effects on this collapsing core. The Rankine-Hugoniot jump condition applied to the case 
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of two colliding cores, indicate that the density in the thin dissipative region where the collision takes place, must 
be increased by about four times the pre-collision density, so we expect an overdensity to be formed at the interface 
region of the colliding cores and also at the core’s central region where gravitational collapse is ongoing. The nature 
of the final distribution of proto-stars will depend on the collision parameters, mainly on the pre-collision velocity 
and the impact factor. 

A word of warning is in order. In any star formation theory, it is very difficult to make numerical calculations with 
realistic initial conditions. Particularly difficult is the case of core collisions, where the parameter space is enormous 
and so we are here forced (like other authors, as already mentioned in the above paragraphs) to use simplified initial 
conditions. Despite this, we assess what are the effects of a collision on the collapsing and colliding cores, as well as 
where and when the accretion centers form, which are the most important issues to be investigated in this paper. 

The outline of this paper is as follows. In Section[2]we describe the initial core, which will be involved in all the 
subsequent collision models. In Section[3]we give the details of the collision geometry and define the models to be 
studied. We shall also describe the modified GADGET2 code and the most important free parameters, which we fix 
by means of a calibration method presented in Section |3.3| In Section [4] we describe the most important features 
of the time evolution of our simulations by means of two-dimensional (2D) and three-dimensional (3D) plots. In 
Sections [5] and [6] we discuss the relevance of our results in view of those reported by previous papers and finally we 
make some concluding remarks. 


2 The initial core 


We consider a uniform density spherical core model with radius Do = 8.0 x 10 1(> cm (= 0.026 pc = 5350 AU) 
and mass Mq = 8 Mq. The average density of this core is po = 9.2 x 10“ ls g cm -3 = 2.3 x 10 6 part cm -3 , 
from which an estimate of its free fall time is £// = ^/37r/(32Gpo) = 6.9 x 10 11 s = 21879 yrs. 

The dynamical fate of a collapsing isolated core is usually characterized by the values of the thermal and rotational 
energy ratios with respect to the gravitational energy, denoted by a and /3, respectively, see | |Bodenheimer et al. (2000)1 
and | Sigalotti & Klapp (2001a) |, | Sigalotti & Klapp (2001b) |. For this paper, all the simulations have been done with 
an initial core having the same initial values a and j3 


a = 


■^therm 

-^grav 


0.3, 


p = 



a) 


As is well known, these values for a and /3 favor the occurrence of core collapse. The speed of sound Co and 
the rotational angular velocity Do have been calculated to satisfy the energy ratios given by Eq. [T] they have the 
following numerical values 


Co = 39956.8 cm s 1 , 

( 2 ) 

Do = 8.09 x 10 13 rads 1 . 

Then the initial velocity of the ith SPH particle is Vi = Do x ry = Do(— yt, Xi, 0). 

Besides, we implement a well known spectrum of density perturbations on the initial particle distribution, so 
that at the end of the simulation they might result in the formation of binary systems. The evolution of a uni¬ 
form density core like the model we consider here has been reproduced by many groups worldwide using different 
codes based both on grids and particles, see for instance |Bodenheimer et al. (2000) |. | Sigalotti & Klapp (2001 a) | and 
| Sigalotti & Klapp (200 1b)| an d the references there in. We have also successfully reproduced the uniform model re- 
sults in |Arreaga et al. (2007)]|Arreaga et al. (2008)| , from which we established the reliability of our calculations for 
the evolution of the core using the GAGDET2 code, as we followed the collapse process until peak densities of the 
order of 2.9 x 10 _1 ° g cm -3 . 

Let us say something about the gas thermodynamics. As the observed star forming regions basically consist of 
molecular hydrogen cores at 10 K with an average density of 1 x 10 -2<) g cm ~ 3 , the ideal equation of state is a good 
approximation. However, once gravity has produced a substantial core contraction, the gas begins to heat. In order to 
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take into account this increase in temperature, we use the barotropic equation of state proposed by |Boss et al. (2000)) . 
Thus, in this paper we carry out all our simulations using the barotropic equation of state 


P = clp 


1 + 



(3) 


where for molecular hydrogen gas the ratio of specific heats 7 is given by 7 = 5/3, because we only take into 
account the translational degrees of freedom of the hydrogen molecule. The most important free parameter is then 
the critical density p cr it, which determines the change of thermodynamic regime. Here we consider only one value 
Of pcrit = 5.0 x 10 -14 g cnr 3 . 

By comparing the results of our previous papers | Arreaga et al. (2007) Arreagaetal. (2008)] with those of 
IWhitehouse & Bate (2006) | for the collapse of an isolated and rigidly rotating core with a uniform density pro¬ 
file, we have concluded that the barotropic equation of state in general behaves quite well and that it captures all the 
essential thermodynamical phases of the collapse, see also [Arreaga & Klapp (2010)| . However, this does not mean 
that the same comparison would be always correct for more complex initial conditions or for collision models such 
as those we consider in this paper. 

Despite these facts and that we know that it is indispensable to include all the detailed physics of the thermal 
transition in order to achieve the correct results, be it in a collision or not, we carry out the present simulations with 
the barotropic approximation, because we know that there are other computational and physical factors that could 
have a stronger influence on the outcome of a simulation. 

The collision cores considered for this work have pre-collision velocities of up to several times the speed of sound, 
see Section [3] and Table [I] For these high velocities a particle agglomeration is formed in the interface between the 
cores and the artificial viscosity used for the calculations transforms kinetic into thermal energy, which increases the 
temperature of the gas. However, as we will show in Section [4~T| the cooling time T CO oi is both much shorter than 
the local sound crossing time, and the free fall time, hence the isothermal condition can be used for the interaction 
region. 


3 Collision models and computational considerations 

3.1 The pre-collision geometry 

The simulations considered in this paper include both head-on and oblique collisions. For the head-on collisions, the 
initial position in cartesian coordinates of the center of mass (CM) of one colliding cores is at (x, y, z) = Xcm 1 = 
( 0 , Ro, Ro) and the other at Xcm 2 = (0, —Ro, Ra )■ Recall that Ro is the initial core’s radius. The oblique collisions 
are characterized by the magnitude and sign of an impact parameter 6 , which depends on Ro. For the 6 > 0 cases, the 
colliding cores are initially displaced along the A'-axis. For the opposite case, when b < 0, the cores are displaced 
again along the A'-axis but in the opposite direction. In all cases the displacement is perpendicular to the symmetry 
axis of the collision (the K-axis). 

For instance, for model M 1 we place the first core so that its center of mass in cartesian coordinates is Xcm ± = 
{Ro/4, Ro, Ro), while the center of mass of the second core is A/m 2 = (—Ro/4, —Ro, Ro)- We consider this 
orientation as a positive impact parameter b, see Tableland Fig. [7] Analogously, for model M3, the center of mass 
of the cores is placed at the coordinates Xcm 1 = (— 6 / 2 , Ro, Ro) and Xcm 2 = ( 6 / 2 , —Ro,Ro), which corresponds 
to a negative impact parameter. 

For the collision models we consider four different impact parameter values 6 = 0, -Ro/4, ± Ro/2, and pre¬ 
collision velocities in the range of 0.23 to 6.4 times the speed of sound. It should be noticed that this range of 
values is well motivated physically, as it follows from a statistical study of core collision in the context of interacting 
galaxies, see [Bekki et al. (2004)} . 

3.2 The pre-collision velocity 

In this paper, we define the pre-collision velocity in the same way for both head-on and oblique collisions, as follows. 
The average velocity of the center of mass of the cores points along the F-axis, and are Vcmx = (0, — V ave , 0) and 
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Vcm 2 = (0, Vave,0)- We define the pre-collision speed V app (with respect to the center of mass of the colliding 
binary system) of the cores by means of V app = 2 V aP e- Therefore, the pre-collision velocity V app of the colliding 
cores, also known as the relative or impact velocity of the cores, should be in the range 0.30 — 20 Mach. 

However, there must be a correlation between the pre-collision velocity of the colliding cores and their separation. 

As we assume here that the cores are initially in close contact, then we can not go to arbitrarily large values of the 
pre-collision velocities. Hence the values used here for V app have been chosen to represent an ample range, where 
the collision timescale can be compared with the free-fall time scale, as discussed in Section[5] 

3.3 The Evolution Code 

We carry out the time evolution of the initial distribution of particles with the fully parallel GAGDET2 code, which is 
described in detail by (Springel (2005)| . GADGET2 is based on the tree-PM method for computing the gravitational 
forces and on the standard SPH method for solving the Euler equations of hydrodynamics. GADGET2 incorporates 
the following standard features: (i) each particle i has its own smoothing length hp, (ii) the particles are also allowed 
to have an individual gravitational softening length e;, whose value is adjusted such that for every time step, ci hi is 
of order unity. GADGET2 fixes the value of Ci for each time-step using the minimum value of the smoothing length 
of all particles, that is, if hmin = min(hi) for i = 1, 2...N, then a = h m in- 

The GADGET2 code has an implementation of a Monaghan-Balsara form for the artificial viscosity, see [Monaghan & Gingold (1983)| 
and |Bals ara (1 995)]. The strength of the viscosity is regulated by setting the parameter a„ = 0.75 and = | x a„, 
see Equation (14) in [Springel (2005)| . We fix the Courant factor to 0.1. The number of smoothing neighbour par¬ 
ticles is fixed to 40 and the allowed variation is fixed to 5 particles, so that the GADGET2 will adjust the kernel’s 
smoothing length such that the number of neighbours is always kept within the range 40 ± 5. 

3.4 The detection of accretion centers 

Let us now describe the modification implemented into the GADGET2 code to detect accretion centers. Any gas 
particle with density higher than p s is a candidate to be an accretion center. We localize all the candidate particles 
for a given time t. We then test the separation between candidate particles: if there is one candidate with no other 
candidate closer than 10 x Tint, then this particle is identified as an accretion center at time t. We define r^t as the 
neighbor radius for an accretion center, given by r^t = 1.5 x hmin . In this way n n t determines a set of particles 
which are within the sphere having this radius and whose center is the accretion center itself. All those particles will 
give their mass and momentum to the accretion center. We change the GADGET2 particle type for all those particles 
within the neighbor radius from being 0 to —1, which is a particle type undefined in GADGET2 and therefore they 
will be not advanced in time anymore. 

We go on to determine an appropriate value of p s by using a uniform and rotating core model studied in previous 
collapse calculations. When we use a low value for p s , for instance, p 3 = 5.0 x 10~ 14 g cm -3 , then too many 
accretion centers are formed even in an early collapse stage, see the right panel of Fig. [2] A better value for p s would 
be p s = 1.0 x 10 -12 g cm -3 , as the evolution of the test model both with and without the implemented code are 
quite similar. 

The dependence of the number of accretion centers on the parameter p s can be considered as a depending way of 
discretize the densest gas region of interest. We emphasize that the radius ri n t and the p s are chosen here such that 
the formation of accretion centers does not affect the subsequent dynamics of the gas outside their influence region, 
otherwise it would be obviously a big problem in our simulations. 

Furthermore, we mention here that [ |Klessen & Burkert (2000)[ have used a clump finding algorithm, fully inte¬ 
grated into the SPH formalism, which makes no use of any parameters like ours p s and r^t- 

However, other methods to detect sink particles like the one described by [Bate et al. (1995)) depends on three 
parameters: (i) a density threshold (which was chosen 10 5 times higher than the average cloud density) to initially 
point to the gas candidate to be a sink particle; (ii) an accretion radius and (iii) an inner radius (which was defined to 
be 10-100 times smaller than the accretion radius). 

[Bate et al. (1995)| applied the sink particle technique to study the collapse of the so called standard isotermal 
test case. Their objective was to extend the collapse calculation beyond the protostar formation by removing all those 
particles in the densest regions, which have very small time steps. 
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[Bate et al. (1995)| applied several collapsing and bounding tests on the gas particles inside the accretion radius 
before they are decided to be accreted or not by the sink particle. Besides, all those particles inside the inner radius 
are accreted by the sink particle regarless of the tests. 

We must clarify that our parameter n n t plays here the role of the inner accretion radius defined by [Bate et al. (1995)) 
and this is why we decide to skip all the tests as well. In a forthcoming revision of our code we will include an exterior 
radius r ex t in order to detect and eventually eliminate all those particles located within an exterior influence region 
of an accretion center. 

[Federrath et al. (2010)| | implemented this sink idea in the mesh based code FLASH, which uses an adaptive mesh 
refiment (AMR) technique. They defined a sink cell and a sink volume by introducing two basic parameters: a density 
threshold and a radius. They went on to apply 6 tests on the sink cell to avoid the formation of spurious sinks. They 
found good agrement between this implementation and that reported by [Bate et al. (1995)) in the case of the core 
collapse. However, they emphasized that the use a sole density threshold parameter for deciding sink creation is 
insufficient when supersonic shocks are present. 

Despite of all this, we carry out the time evolution of all the initial distribution of particles for the colliding models 
with a fixed value of p s = 5.0 x 10~ 14 g cnF 3 , because of the failure to evolve many of the particles of the colliding 
models until densities higher than p s , otherwise the time step of the GADGET2 code becomes very small when not 
many particles are able to reach densities beyond p s . As expected, all the collision models show a strong tendency 
to collapse like the isolated core, as can be seen in Fig. [3] The application of the limited technique described in this 
Section makes sense for us, as the main aim of this paper is only to detect the places where the densest region are 
formed when two rotating cores collide with a moderate pre-collision velocity. 

3.5 Resolution 

According to [Truelove et al. (1997)} , the resolution requirement for avoiding the growth of numerical instabilities is 
expressed in terms of the Jeans wavelength A j, which is given by 


A j = 



(4) 


where G is Newton’s universal gravitation constant, c is the instantaneous sound speed and p is the local density. To 
obtain a more useful form for a particle based code, the Jeans wavelength A j is written in terms of the spherical Jeans 
mass Mj, which is defined by 



Nowadays it is well known that the Jeans requirement l < Aj/4 (where l is a characteristic length scale of the 
grid for a mesh based code) is a necessary condition to avoid the occurrence of artificial fragmentation. For particle 
based codes, [Bate & Burkert (1997)) showed that there is also a mass limit resolution criterion which needs to be 
fulfilled besides that of [Truelove et al. (1997)| . They showed that an SPH code produces correct results involving 
self-gravity as long as the minimum resolvable mass is always less than the Jeans mass Mj. 

Now, if we approximate the instantaneous sound speed by c = s/p/ P, then according to Eq. [ 3 ] we have that 

3 
2 

( 6 ) 

Following [Bate & Burkert (1997)| , the smallest mass that an SPH calculation can resolve is m r ~ Mj/(2N ne i g h), 
where N ne i g h is the number of neighbor particles included in the SPH kernel. For our collision models to comply 
with the Jeans requirement, the particle mass m p must be such that m p /m r < 1. 

In this work, we set N = 1, 000, 000 SPH particles for representing the initial core configuration in each model. 
By means of a rectangular mesh we make the partition of the simulation volume in small elements each with a volume 
A* Ay A z\ at the center of each volume we place a particle -the i — th, say- with a mass determined by its location 


Mj = 
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according to the density profile being considered, that is: rm = p(xi,yi, Zi ) Air Ay A z with i = 1, N. Next, we 
displace each particle from its location a distance of the order Ar/4.0 in a random spatial direction. 

For verifying that the Jeans stability condition is satisfied we will use the collision model HO — 3, which is the 
one that reached the highest maximum density of all models, see Table[7] For this model we follow its collapse until 
a peak density of pmax = 2.8 x 10 -ll] g cm - ' 5 . The particle mass is m p = 2Mq/N p = 8 x 10 -6 Mq, where 
Mt — 16Mq is the total mass contained in the simulation box and N p is the total number of particles, which are 2 
million for all the models in this paper. 

The minimum Jeans mass for model HO — 3 is given by (Mj) HO _ 3 ~ 8.6 x 10~ 4 M@, from which we obtain 
m r = 1.0 x 10 -5 Mg. Thus, for model HO — 3 we obtain the ratio m p /m r = 0.74, and the Jeans resolution 
requirement is satisfied very easily. 

For the rest of the collision models, the minimum Jeans mass is expected to be greater than for model HO — 3 
because their maximum density is less than for model HO — 3. It is then clear that for all models the Jeans length 
criterion is satisfied as well. 


3.6 Collision models 

We summarize the collision models considered in this paper in Table |T] whose entries are as follows. Column one 
shows the label chosen to identify the model. Column two indicates the impact parameter b of the collision models 
in terms of Ro', note the appearance of a sign before the magnitude of b, which is related to the orientation of the 
colliding cores along the X-axis, see Section [3T] Column three indicates the pre-collision velocity of the colliding 
cores, expressed in terms of the initial speed of sound in the core, whose numerical value is given in Eq. [2] In 
the fourth column we show the number N acc of accretion centers found for each system, both for the case of an 
isolated core and for the colliding models. In the fifth column we show the average number of SPH particles Np acc 
per accretion center formed in each model while in the sixth column we show average mass of the accretion center 
Mav/M aurl . As a way of comparing our simulations with other simulations, in the seventh column we show the peak 
density pmax reached in each run. 


4 Results 

In order to show the results of our simulations, we present iso-density plots for a slice of matter parallel to the X — Y 
plane. The procedure is as follows. We first locate the SPH particle with the maximum density in the entire volume 
space of a simulation. The z-coordinate of that particle, say z ma x, determines the height of the thin slice of material. 
The width of the slice is determined so that about 10, 000 SPH particles enter into the slice, which is centered around 
the Zmax coordinate. 

Once the SPH particles defining the slice have been selected, we set a color scale related to the iso-density curves: 
yellow indicates areas with higher densities, blue indicates those with lower densities, and green and orange indicate 
intermediate density regions. It should be noted that there is no relation between the density colors associated with 
different panels even in the same plot. At the bottom of each iso-density panel, we include two numerical values to 
illustrate the different stages of the evolution process: the left value is the peak density p(t) at time t; and the right 
one is the actual time t in seconds. 


4.1 The general picture of a core collision 

Let us now briefly describe the general picture of a collision system by using the evolution of model HO — 1 as an 
example, calculated using the normal GADGET2 code just for comparison. In the first top left panel of Fig. [4] one 
can see that two identical cores are just placed one against the other at t = 0. Each core shares those particles which 
are closer to the edge of the neighboring core. 

As each core in the colliding system is collapsing separately as well, we see in the right bottom panel of Fig. [4] 
the formation of a filament in each central core region. Such filaments are indeed connected by a gas bridge that is 
formed in the interface region of the cores. 
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The interface particles are compressed and a front perturbation is formed that propagates into the two cores. The 
artificial viscosity then transforms kinetic energy into heat. We expect that the barotropic equation of state proposed 
by |Boss et al. (2000)] and described in Eq.[3]mimics the radiative cooling needed as this heat must be radiated away 
in a very short timescale, so that we can assume that the cores and the interface region between the cores remain 
isothermal. The perturbation wave that propagates into the cores can disturb the filaments in the central region. 

During the collision process a slab of material is formed along the collision front. The cores then expand but 
eventually collapse again. The compressed slab formed is susceptible of having various instabilities, which have been 
studied by (Vishniac (1983)| for the linear regime, and by [ [Vishniac (1994)| for the non-linear case. For the present 
work, the relevant instabilities are the shearing and the gravitational. The shearing instability dominates at low density 
while the gravitational one takes over for densities above the critical density 

PsHear ~ 7.86132 X lO " 22 (j^ (J) 0) 

where i ) is the amplitude of the perturbation, A ntsi the length of the unstable mode, and T the temperature, see 

(Heitsch et al. (2008)| . 

Of the shearing instabilities, the main ones are the non-linear thin shell instability (NTSI), and the Kelvin- 
Helmholtz (KH) instability. The NTSI instability is expected to occur in the shocked slab just after the core collision, 
and the non-linear bending and breathing modes could also be present. From our numerical results we estimate that 
the parameters in Eq.[7]are A ntsi ~ 0.1 pc, 77 ~ 2 A ntsi and T ~ 10 K, and so the NTSI and KH instability 

are suppressed when the density goes above p cr it ~ 1.57 x 10~ ls g cm -3 . 

In the next sections we shall describe the results obtained with the modified GADGET2 code for the time evolution 
of all models. The aim is to locate accretion centers and to determine their properties along the filaments and at the 
interfacial region as well. 

4.2 Head-on collisions 

In Fig. [5]we show the last snapshot available for each collision model, or in other words, for the sake of brevity, we 
omit presenting the earlier states of the evolution, like the first three panels shown in Fig. [4] 

With V app slightly less than Co, which corresponds to model HO — 3, the initial cores are slightly disturbed, as 
the central filament shows a little stretching. 

Most of the gas in the filament remains there as the collapse of each core progresses further. As we see in Fig. [2] 
many accretion centers are expected to be formed along these filaments, giving place to the formation of a peculiar 
chain of protostars when the gas cools further. 

In the top panel of Fig. [ 6 ]we show the velocity field of the particles located near the filament. It can be seen that 
these particles flow into the filament. In Fig. [ 8 ] we show a 3 D visualization of the same filament of model HO 2 seen 
from a different angle; in this view, it is possible to observe that a very small over-density is formed at each filament’s 
end, at the colliding interface between the two collapsing cores. The origin of these new over-densities is the density 
perturbation caused by the core collision. 

With a higher pre-collision velocity, self-gravity will have less time to act on the core and many more particles 
from the cores would flow very quickly to the cores’ interface, up to the point that the two initial cores can no longer 
be distinguished, as can be seen in models HO — 6 and HO — 7. For these models, the net effect of the collision 
is then to replace the two initial colliding cores by only one core, which rapidly forms a long central filament-like 
structure. In the bottom panel of Fig. [ 6 ] we observe again the formation of a strong filament in the interface of the 
cores that expands rapidly as a strong flow of particles go by the ends along the X-axis. However, the expansion 
eventually stops and the filament bounces to re-initiate the expansion in the transversal direction, along the Y -axis. 
This behavior has previously been observed by (Anathpindika (2009a)| . 

Thus, one effect of having a higher pre-collision velocity is that the density perturbation is stronger as can be seen 
in the last panel of Fig. [5] If V app is high enough, we see that the net effect of the collision is to replace the two initial 
colliding cores by one core which may not be collapsing anymore. 
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4.3 Oblique collisions 

The outcome for the oblique models is illustrated in Fig. [7] where we show only the last available snapshot of each 
model. As expected, the results for model Ml, with a low pre-collision velocity and with small impact parameter 
value, are quite similar to those obtained for the head-on collision models HO — 1, HO — 2 and HO — 3. Hence, 
it is unnecessary to calculate more models with increasing V app with a small impact parameter, because we would 
expect to see similar results to those obtained for the head-on collision models. 

We then go on to models M2,M'21 and M 22, where the impact parameter has been increased to 6 = Rq /2, again 
with the positive orientation. The pre-collision speed is low for the two former models and higher for the later model. 
Similar results are obtained for models M 3 and A/4, where we change the orientation of the impact parameter. 

The net effect of a high V app on the colliding system is evident in models M22 and A/4, as the two cores get 
dispersed to give place to a long filament. It must be emphasized that the orientation and the width of these filaments 
are different from the ones obtained for the high V app head-on collision models. 

4.4 Formation and location of the accretion centers 

In order to show the results of this section, we present 3 D plots of all the particles located in a central cube (with half 
the length of the side of the original simulation box) and whose density is greater than 1.2 x 10 1 ' g cm -3 , see 
Fig. [ 8 ] As in Section [4] we again set a color scale related to density: red and orange indicate particles with higher 
densities, blue those with lower densities, and and yellow the intermediate density particles, see the caption of Fig. [ 8 ] 
It should be noted that there is no relation between the density colors associated with different panels even in the same 
plot. We rotate each panel as needed for the purpose of better appreciating the location of the accretion centers. 

As expected, the pre-collision velocity V app also has a strong influence on the number and location of accretion 
centers, as it seems that the greater the V app the lower the N acc , As can be seen in Fig. [ 8 ] if V app is too low, the initial 
binary core system survives even to an advanced stage of the evolution, and the accretion centers are mainly localized 
in the core’s central filament. But if V app is too high, then the initial binary core system is destroyed by the strong 
flow of particles, giving less chance for accretion center to form. For models HO — 3 and HO — 6 , the accretion 
centers are indeed located mainly in the cores’ interface, as can be seen in the close-up shown in Fig. [ 8 ] 

It is therefore remarkable to realize that models HO — 2 and HO — 3 have accretion centers forming outside the 
core’s central filament, as can be seen in Fig. [ 8 ] As we can see, there must be a small range of velocities V app which 
favors the formation of accretion centers in the gas bridge connecting the two cores, see Fig. [9] 

From Table [T] we first notice that the oblique collision models have a larger number of accretion centers than the 
head-on models. It seems that more accretion centers are created if self-gravity has more time available to act on the 
collapsing cores. 


5 Discussion 

For this paper we carried out a fully 3 D set of numerical hydrodynamical simulations within the framework of the 
SPH technique, aimed at following the formation of accretion centers in a collision process of two small rotating 
cores with uniform density. 

The initial conditions for the isolated core are such that it will collapse because its initial (a + /3) CO re < 0.5. 
However, when we consider the translational kinetic energy that comes from the pre-collision velocity V app of the 
cores, the collision system is expected to be slightly unbound, that is, (a + /3) C0 iu s i 0 n > 0.5 for all models. As a 
result of the collision process, a large fraction of the translational kinetic energy will be transformed into heat, which 
will be radiated away because the cores and the perturbation front that forms in the interphase between the cores are 
isothermal. This must somehow produce a reduction of the global value of (a + /?) for our collision models, as we 
observe that most of them are on the verge of gravitational collapse, see Fig. [3] 

We have seen that the result of a collision simulation depends mainly on two physical factors that compete to exert 
influence on the collision process. These factors are the self-gravity and the flow of the particles from the colliding 
interface. The occurrence and extension of these factors are regulated by the value of the pre-collision velocity V app 
of the collision system. 
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With a very slow pre-collision velocity, self-gravity dominates the evolution of both cores. In such cases, the main 
effect of self-gravity on the system is to reduce the core size. The collision system quickly reaches a configuration 
in an advanced state of gravitational collapse, consisting of the remains of each collapsed core linked by a bridge 
of matter. We note that the accretion centers in the simulations are always located in the gas remnants of the cores. 
In most collision cases, we expect these accretion centers to separate from the filamentary structure, so that each 
fragment is expected to evolve towards virial equilibrium. 

On the other hand, when a collision takes place with a significant pre-collision velocity, self-gravity plays a minor 
role and the fate of the system is entirely determined by the strong flow of particles in the interface of the cores. In 
these cases, the main effect on the system is to replace the two initial collapsing cores by a single core that rapidly 
gets perturbed to give place to an elongated and irregular filamentary structure formed in the interaction region. This 
filamentary structure increases its size very rapidly as it is continuously fed by in-falling material coming from the 
original cores. In these cases, we observe that accretion centers are formed in the central region of this filament 
structure. 

However, in view of the limited nature of our simulations, the future of this filamentary structure is as yet unclear; 
but it is likely that most of the gas already in the filament will stay in this filamentary structure as the collapse 
progresses further. Although accretion centers may be formed near the bridge edges, as we see in the case of model 
HO — 2, these accretion centers will be very thin and small, so it could be the case that their diameter will be greater in 
general than the Jeans length. This is another reason why the bridge would not be an appropriate site for proto-stellar 
formation until it cools further. 

The results obtained for model HO — 3 have been emphasized by zooming in on the filament’s central region, as 
shown in Fig. [9] We here see the formation of a well defined central accretion center, which is surrounded by many 
smaller accretion centers. The properties of this clustered distribution of protostars will be discussed elsewhere. 

Now, according to Table[2] the maximum evolution time reached by the models is in the range 0.6 < t evo i/tff < 
1.3, where the free fall time f// was defined in the first paragraph of Section[2] Except for models HO — 6 and 
HO — 7, all the other models have reached a t evo i greater than tff. It should be noticed that most models have 
therefore evolved in time up to an advanced stage of collapse. Unfortunately, it is not possible to compute most of the 
models any further, as the time step becomes very small. 

The relevant time scale for the collision is given approximately by t co i ta Ro/V app . For the models of this paper, 
this collision time is in the range 0.5 < t co i/tff < 12. For model HO — 6 these time scales are of the same order, 
tcoi/tff ~ 1. By looking at Fig. [5] we see that the pre-collision velocity of model HO — 6 is the maximum collision 
velocity that allows an initial binary colliding core system to remain as a binary system capable of undergoing further 
gravitational collapse and eventually of producing multiple protostars. It is in this sense that our results have to be 
compared with those of [ [Whithworth & Pongracic (1991)} , who found that it is more likely that a colliding system 
will in general result in disruption and dispersal of the core involved, with no chance of forming protostars. 


6 Concluding Remarks 

In this paper we have been able to demostrate some of the essential features of the two core collision processes. In 
a first approximation, we have investigated the number and location of the accretion centers formed in this collision 
scheme. 

The modified GADGET2 code tested here and described in Section[33]is a first step towards the full sink-particle 
technique implementation already in use worldwide. Our code is still incomplete as we have not yet included all the 
particle tests suggested by |Bate et al. (1995)) and |Federrath et al. (201 0) |. However, as we have applied the same 
code to all the models, a comparison is possible. The main results of this comparison are shown in Table [I] and in 
Fig. [8] 

It must be taken into account that (i) we have used a low density threshold p s for a particle to become an accretion 
center and (ii) there is no restriction on the minimum number of particles that must be captured for a particle to 
become an accretion center: so we have seen that there are a large number of accretion centers having just very few 
captured particles. For these two reasons we must expect an overproduction of accretion centers in our simulations. 
Finally, we recall that accretion centers are gas particles and can collide or merge. In view of this, the future of the 
accretion centers is still unclear, and requires taking the calculations further in time. 
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Despite this, it is interesting to note that not all the accretion centers in the simulations were always found located 
in the remnants of the cores. There were also accretion centers found in the bridge of matter, as was the case for 
model HO — 2; and in the central collision region, as occurred for models HO — 3 and HO — 6. Hopefully all these 
accretion centers will end up as proto-stellar seeds forming either a peculiar chained structure or a clustered structure, 
respectively. For these reasons, our hope is that the results presented here can be considered as further evidence that 
core collisions may have an important influence on the star formation process. 
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Figure 1: The initial geometry of the colliding cores for head-on collisions (left panel), and for oblique collisions with 
a positive impact parameter b (right panel). 



Figure 2: The last snapshots available for the collapse of the isolated core, for which we use the modified GADGET2 
with density cutoff given by p s = 1.0 x 10 -12 g cm -3 (left panel), and p s = 5.0 x 10 -14 g cm~ 3 (right panel). 


Table 1: The collision models. 


Model 

b 

Ro 

Vgpp 

Co 

Nacc 

N Pace 

M av /M Q 

Pmax [gr/cm 3 ] 

MCmR 



27 

52 

5.80 x 10" 5 

1.17 x 10 1U 

HO-1 

0 

0.23 

80 

22 

1.97 x 10" 5 

7.05 x 10 -12 

HO-2 

0 

0.46 

29 

5 

4.20 x 10" e 

6.7 x 10- 11 

HO-3 

0 

0.92 

30 

6 

5.0 x 10" y 

2.8 x 10~ 1U 

HO-5 

0 

3.9 

104 

5 

3.9 x 10" H 

3.37 x 10 ii 

HO-6 

0 

2.9 

37 

5 

4.0 x 10" B 

1.0 x 10~ 1U 

HO-7 

0 

5.7 

9 

57 

4.6 x 10" 5 

1.66 x 10" 13 

Ml 

1/4 

0.23 

210 

7 

5.9 x 10" y 

1.2 x KT 11 

M2 

1/2 

0.25 

103 

15 

1.2 x 10" 5 

7.0 x 10” 12 

M22 

1/2 

6.4 

42 

5 

4.4 x 10" a 

7.1 x 10” 11 

M21 

1/2 

0.51 

214 
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4.2 x 10" 6 

1.8 x 10” 11 

M3 

-1/2 

0.51 

156 

5 

4.3 x 10" y 

1.8 x 10” 11 

M4 
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6.4 

30 
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4.0 x 10 -6 

1.2 x KT 11 
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Figure 3: Time evolution of the peak density for: the isolated core (top panel), the head-on collision models (middle 
panel), and the oblique collision models (bottom panel). 
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Figure 4: Iso-density plot illustrating the gravitational collapse process for model HO 1 computed with the normal 
GADGET2 code. 
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Figure 5: Isodensity plots for the last available snapshot in the head-on collision models HO — 1 (top left panel), 
HO — 2 (top middle panel), HO — 3 (top right panel), HO — 6 (bottom left panel), HO — 5 (bottom middle panel), 
and HO — 7 (bottom right panel). 
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Figure 6: Velocity field of the filament in models H 02 (top) and HOG (bottom), see the top middle panel and bottom 
left panel of Fig. [5] respectively. 
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Figure 7: Iso-density plots for the last available snapshot for the oblique collision models Ml (top left panel), M2 
(top middle panel), M 21 (top right panel), M 22 (bottom left panel), M 3 (bottom middle panel), and M4 (bottom 
right panel). 
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Figure 8: Location of the accretion centers for the last available snapshot of models HOI (first line left), HO 2 (first 
line right), HOS (second line left), H05 (second line right), H06 (third line left), M22 (third line right), M3 (fourth 
line left), and M4 (fourth line right). The colors are a density scale for log 10 according with the following 
ranges: blue ( < 2), green (2 — 3), yellow (3 — 4), orange (4 — 5) and red (> 5). 
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Figure 9: Zoom in of the left panel of the second line in Fig. [8] where the accretion centers are shown for model 
H03. 


Table 2: The evolution and collision time of the collision models. 


Model 

tevol/t f f 

tcol/tff 

HO-1 

1.3 

12.6 

HO-2 

1.3 

6.3 

HO-3 

1.1 

3.15 

HO-5 

0.82 

0.74 

HO-6 

0.9 

1.0 

HO-7 

0.75 

0.5 

Ml 

1.3 

12.6 

M2 

1.3 

12.6 

M21 

1.3 

5.6 

M22 

0.74 

0.45 

M3 

1.3 

5.6 

M4 

0.92 

0.45 
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